function diff = rep_cost_p_fun( m, param)

    tmp2 = -2.*param.c.*(1-param.alpha)./param.sigma^2./(param.gamma_R_2-param.gamma_R_1).* ...
        ((1./(1-param.alpha)-param.gamma_R_1).*param.m_u.^param.gamma_R_1.*int_cost1_fun(param.m_u,param).*param.gamma_R_1./param.m_u ...
        -(1./(1-param.alpha)-param.gamma_R_2).*param.m_u.^param.gamma_R_2.*int_cost2_fun(param.m_u,param).*param.gamma_R_2./param.m_u ...    
        +(param.gamma_R_2-param.gamma_R_1).*(param.zeta+delta_fun(param.m_u,param))./param.m_u );

    tmp1=J_0_fun( param.m_h, param);

    A = (tmp1.*(param.m_u./param.m_h).^param.gamma_R_2.*param.gamma_R_2./param.m_u - tmp2) ./ ...
        ((param.m_u./param.m_h).^param.gamma_R_1.*param.gamma_R_1./param.m_u - (param.m_u./param.m_h).^param.gamma_R_2.*param.gamma_R_2./param.m_u);

    diff = -2.*param.c.*(1-param.alpha)./param.sigma^2./(param.gamma_R_2-param.gamma_R_1).* ...
        ((1./(1-param.alpha)-param.gamma_R_1).*m.^param.gamma_R_1.*int_cost1_fun( m,param).*param.gamma_R_1./m ...
        -(1./(1-param.alpha)-param.gamma_R_2).*m.^param.gamma_R_2.*int_cost2_fun( m,param).*param.gamma_R_2./m ...    
        +(param.gamma_R_2-param.gamma_R_1).*(param.zeta+delta_fun(m,param))./m ) ...
        + A.*(m./param.m_h).^param.gamma_R_1.*param.gamma_R_1./m ...
        - (tmp1+A).*(m./param.m_h).^param.gamma_R_2.*param.gamma_R_2./m;

    diff = - diff.*(m>=param.m_h);

end